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In earlier work P, @, [3| a systematic extension of the van der Waals square gradient model to 
non-equilibrium one-component systems was given. In this work the focus was on heat and mass 
transfer through the liquid-vapor interface as caused by a temperature difference or an over or 
under pressure. It was established that the surface as described using Gibbs excess densities was 
' in local equilibrium. Heat and mass transfer coefficients were evaluated. In our first paper Q we 

' discussed the equilibrium properties of a multi-component system following the same procedure. 

, In particular, we derived an explicit expression for the pressure tensor and discussed the validity 

' of the Gibbs relation in the interfacial region. In this paper we will give an extension of this 

^ ' approach to multi-component non-equilibrium systems in the systematic context of non-equilibrium 

Q , thermodynamics. The two-dimensional isotropy of the interface is discussed. Furthermore we give 

h-^ • numerically obtained profiles of the concentration, the mole fraction and the temperature, which 

' illustrate the solution for some special cases. 

i> : 

I. INTRODUCTION. 

o . 

c/2 Because of the lack of local equilibrium the extension of non-equilibrium thermodynamics to a continuous description 

' of an interface is not straightforward. In earlier work HI, [3, II , coauthoredby one of us (DB) , we were able to show that 
' such an extension was possible for one-component fluids, with all the variables dependent on the normal coordinate 
^ for a planar interface. Temperature gradients, pressure differences and the resulting heat flux and evaporation or 
condensation fluxes were determined through the interface. For systems away from equilibrium, square gradient 
models have been used before, we refer to [1,0, 0| in this context. Very little work has been done on systems with a 
varying temperature 0, and on two-phase systems, however. The systematic treatment of heat and mass transport 
Q I through the liquid- vapor interface, along the lines sketched in [H, 0, 3, was to our knowledge, new. In the first paper 
Q we extended the analysis to equilibrium multi-component systems. Explicit expression for all the thermodynamic 
I I quantities were given. In particular, we derived the expression for the pressure tensor and discussed the validity of the 
^ . Gibbs relation in the interfacial region. In this paper we will extend this approach to multi-component non-equilibrium 
' mixtures in three-dimensions. Because of considering three-dimensional systems it is possible to derive systematically 
, all system properties, which were postulated in the one-dimensional description. 

The traditional local equilibrium hypothesis implies that in non-equilibrium a small volume of the system at each 
moment of time can be considered to be in equilibrium. Thus, all thermodynamic relations, valid for the whole system, 
remain valid for this small volume. The important assumption is that the state of each small volume is determined 
only by the properties of this volume, and no other ones. To describe the interface in equilibrium, one needs to 
. introduce a dependence on the density gradients however. Such a description is not local in the traditional sense: 
the system behavior in a small volume depends on the properties of this volume and the properties of the nearest 
neighborhood. This implies that one cannot apply the usual local equilibrium hypothesis to the interfacial region. 
We will follow the traditional procedure used in non-equilibrium thermodynamics. In Sec. [11] we review briefly the 
^ , main results of the equilibrium square gradient model for the interface. We extend all equilibrium results, derived 
' in [J to non-equilibrium. Explicit expressions for the thermodynamic properties of the non-equilibrium mixture are 
given. To be able to derive the entropy production one needs the Gibbs relation. In was shown in 4] that within 
the interfacial region one can speak about two kinds of the Gibbs relation, the ordinary Gibbs relation and the 
spatial Gibbs relation. We show in Sec. [IT] how these two static relations are combined to the non-equilibrium Gibbs 
relation. For a one-component system the non-equilibrium Gibbs relation reduces to the one, given in [l|. Within 
the framework of the one-fluid approach we discuss the expression for the homogeneous Helmholtz energy /o which 
is needed to determine thermodynamic quantities in Sec. [Illj . In Sec. jlVj we give the hydrodynamic equations. We 
use the so-called one-fluid approach which has been shown [M, to be appropriate for common mixtures. In Sec. [V] 
we discuss the consequences of the special surface symmetry. After deriving the entropy production and using the 
Curie principle, we give the linear relations between the thermodynamic fluxes and forces in Sec. |VIj . It is found 
that for instance the resistances for transport through and into the interfacial region will in general contain square 
gradient contributions. In Sec. [Vllj we give for some examples numerically obtained concentrations, mole fraction 
and temperature proflles for a binary mixture with stationary mass and heat transport through the interfacial region. 
Concluding remarks and a discussion are given in the last section. 
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II. THE SQUARE GRADIENT MODEL. 



A. Equilibrium surface. 

A mixture can be described by mass densities and temperature T as independent variables. We will use p(r) as 
the total mass density of the mixture, or t;™(r) = l/p(r) as the mass specific volume, and {^™(r), . . . , ^™_]^(r)} as 
the mass fractions of components. Furthermore we will write ^"^ instead of the set of arguments {'^™, . . . , C^i} ^-nd 
V^™ instead of the set {V^™, . . . , y£J^_i} to simplify the formulas. 

Assuming that the specific Helmholtz energy can be written as 

/™(r) = /^-(T, p, r ) + /c'"(p, r , vp, vr ) (n.i) 

where 

/C"(p, r, Vp, vr) ^ 2^ + ^ ^ Vp(r).Vr (r) + 2 E ^ VC(r)-Ve™(r) (II.2) 

and where all coefficients k, and Kij are assumed to be independent of the temperature, one can derive the chemical 
potential /i„ of the n'^ component and the reduced chemical potential V'fc = Mfc — Mm where k = l,rt — 1 (all integers 
from 1 to n — 1), which are constant through the surface in equilibrium. We refer to ^ for the expressions for the 
chemical potentials, p(r) and other quantities. These quantities are related by the ordinary thermodynamic relation 

r(r) = /i™ - p (r) ^;™(r) + ^ ^fC™ (0 (11.3) 

i=i 

Other thermodynamic potentials, like the specific internal energy 

M'"(r) = /™(r) + s^Cr) T (II.4) 

The entropy and the enthalpy are also given. 

In the interfacial region, the pressure has a tensorial behavior: 



o-Q/3(r) = p (r) 5al3 + 7a/3(l") (II. 5) 



where the tension tensor is given by 



,.^^.^ dp{v)dp{v) ^ / 94'"(r)gp(r) ^ dp(v)d^T{^) \ , y\.^,n ^^(r) ^Cr(r) 

"'^ dxR ^ * V cJxq dxR dxa dxg J ^ dxa dxg 

^ z— 1 ^ ^ — 1 

For a flat surface pir) is the parallel pressure. 

An important part of the equilibrium description is the relation between the rate of change of thermodynamic 
potentials and the independent thermodynamic variables, in other words the Gibbs relation. In the interfacial region 
properties may change with the position, so we have to speak about two kinds of Gibbs relations: ordinary Gibbs 
relations and spatial Gibbs relations. Ordinary Gibbs relations relate the changes of thermodynamic variables at the 
given point in space for different states. For the internal energy the ordinary Gibbs relation was found to be 

n-l 

<5w™(s"(r),«™(r),r(r)) - T(5s"(r) + ^ ^"^^(r) ~p{v)5v^{v) (II.7) 

i=l 

Spatial Gibbs relations relate differences of thermodynamic variables for the given state at neighboring points in 
space. For the internal energy it was found to be 

VM™(r) =rVs™(r)+^Vz'"VC(i-)-p(r)V?;™(r)-v"(r)^^^ (II.8) 



Notice in particular the last contribution on the right hand side. 
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B. Non-equilibrium surface. 



In order to describe non-equilibrium processes in a multi-phase mixture using thermodynamics, one must assume, 
that all thermodynamic quantities are defined in each point in space and at all times. This in particular is also in the 
interfacial region. In non-equilibrium the density p(r, t) and mass fractions ^'"(r, t) depend on the time explicitly. 
We do not have the restriction of a constant temperature and chemical potentials: r(r, t), /x™(r, t) and ij^Y^ir, t) 
may therefore depend both on position and time. New variables which appear in the non-equilibrium description are 
the velocities of each component {vi(r, t), . . . , v„(r, t)}. We will however use the barycentric velocity of the whole 
mixture 

1 " 

v(r, t) = ^^(77) E ^)^^ ("-9) 

and the diffusion fluxes 



j;"^per(vfc-v), fc = i,n-i (11.10) 

as independent variables. 

To describe non-equilibrium inhomogeneous systems we shall assume that all the relations between thermodynamic 
variables valid in equilibrium, which were discussed in the previous section and in more detail in , remain valid away 
from equilibrium. Away from the surface this is the usual assumption made in non-equilibrium thermodynamics. In 
the interfacial regions it extends this assumption to places where the gradient contributions become important. As is 
said in 10], the validity of such a hypothesis can be verified only by experiment. 

An important part of the local equilibrium hypothesis is the relation between the rate of change of thermodynamic 
potentials and independent thermodynamic variables, in other words the Gibbs relation. Similar to the description 
of a homogeneous fluid, we extend the equilibrium Gibbs relations in the simplest way. One needs to make one 
important observation before such an extension, however: equilibrium equation Eq. (|II.7|) describe the change of local 
thermodynamic variables between two different states at a fixed point in space. These two states can be separated in 
time. So we can say that this equation describes the change of local thermodynamic variables in time at a fixed point 
in space: 

nr. „ ^ ™ - 'i; .) ^ , ^ „,„) 

i—l 

We similarly use the equilibrium spatial Gibbs relation for the specific internal energy, Eq. pi.8[) . for non-equilibrium 
case: 

T(r, t) Vs"(r, t) = Vu"(r, t) - i^Ti^, t) VC(r, t)+p{v, t) Vt;™(r, t) - v^{v, t) ^22i^ElA (11.12) 

See also relevant discussion in Appendix [X]. Contracting Eq. (|II.12[) with v(r, t) and summing with Eq. (jlLlip we 
find 

n-1 



- -^^ir^ 5: V^r(r, t) + pir, t) --^^ - v-ir, t) v,(r, t) ^il^ (11.13) 

i—l ^ 

where we used the substantial (barycentric) time derivative 

Eq. (jII.13|) is the Gibbs relation for the non-equilibrium two-phase mixture including the interface. One can show 
that it reduces to the Gibbs relation used by Bedeaux et al 1] for the case of a one-component fluid. The above 
analysis gives more insight in the origin of the contribution proportional to the divergence of the surface tension field. 
This was not clarified in the analysis of Bedeaux et al. This defined all quantities and gave all relations needed in the 
non-equilibrium description using a generalization of the hypothesis of local equilibrium. We will further omit the 
arguments (r, t) to simplify the notation. 
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C. Mixing rules for the gradient coefficients. 

Molar coefficients k™, and are related to the volume coefficients n^j as following (see Q): 

n—1 n—1 

^ y: cr Kj + 2 e + k" 

n-l 

C =P E erk^.+pk^ (11.15) 

_ 2 T^-u 



where 

(11.16) 

All these coefficients are in principle known functions of the densities. In practice they are not known for mixtures. 
Only the values for pure components are more or less known. Thus, it is necessary to express the cross-coefficients 
in a form such that we can approximate them using pure-component values. From the equilibrium analysis for the 
specific quantities per unit of volume [3| one can see, that are simply related to the pure-component coefficients. 
We will assume then that, for instance, = is the coefficient for the pure component i. Cross-coefficients can 
then be approximated by one of the so-called mixing rule for the gradient coefficients. We will assume the following 
mixing rule 

<3 = ("-17) 

This mixing rule for the volume coefficients is analogous to the mixing rule for given in Sec. . 



■\.v — V 

k = 



III. HOMOGENEOUS HELMHOLTZ ENERGY OF ONE-FLUID MIXTURE. 

To obtain the homogeneous specific Helmholtz energy we use the most common one-fluid equations of state: 

- van der Waals equation of state: 

r-r MT, g) RTc Airj. 2 nui\ 

V — B[£^) l — B{Qc 

- Soave-Redlich-Kwong equation of state: 

(T A{T,0 _ RTc A{T,0 ^ 

Ps.k[^.c,i;) viv-B{0) 1-BiOc 1-B(0c'' ^ ^ 

- Peng-Robinson equation of state: 

(rr. ^ _RT A{T, g) ^ RTc A{T, 2 (.rr n^ 

^^"^ ' v-B{0 v^ + 2B{0v-B^0 l-S(C)c l + 2S(e)c-B2(^)c2 ^ ' 

where ^ is a short notation for the molar fractions {^1, . . . , S,n}- 

In the one-fluid approach, constants A{T) and B, depend on the fractions of the species due to the mixing rules: 

n 

A(r, a, ...,?«) = E a,k{T)^^^k 

''n=^ (III.4) 

S(a, ...,^n) = E ^fe^ 

fc=l 



where usually aik{T) — ^ ai(T) ak{T), and ak{T) and bk are the corresponding coefficients for pure substances. The 
matrix aik(T) is symmetric in it's indexes. 
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To find the homogeneous specific Helmholtz energy fo(T, c, ^) we integrate the equations of state over the volume 
at constant temperature and molar fractions of the components. The integration constant should be chosen such 
that the specific Helmholtz energy of the system with a small mixture's concentration (c ^ 0) is equal to the specific 
Helmholtz energy for a mixture of ideal gases. Integrating the equation of state and using Eq. (jC.7|) in Appendix [C] . 
we obtain following expression for the homogeneous specific Helmholtz energy: 

/o(T, c, = -RT In (-|-^^:^ (1 - B{0 c)) - ^(T, c^(B(e) c) (III.5) 

where (^(w) has corresponding expression for each equation of state: 

'PsRKi^) = (III.6) 

IV. HYDRODYNAMICS OF ONE-FLUID MIXTURE. 

We can now derive all hydrodynamic equations, using the conservation laws of matter, momentum and energy. The 
laws of conservation of mass can be written as 

(iv.i) 

=-V-Jr, k=l,n~l 
The momentum conservation law, or the equation of motion can be written as 

where g is the gravitational acceleration, aap is the thermodynamic pressure tensor defined by Eq. pi.5[) and 7rQ,^(r, t) 
is the viscous pressure tensor, which still is to be determined. The viscous pressure tensor without subscripts will be 
indicated by H(r, t). We assume, that this tensor is symmetric. 
The law of energy conservation is (see 

r/p™ 

p +V(Je-pve") = (IV.3) 
at 

where the total specific energy e™ is given by 

e'"(r, t) = w"(r, t) + T'"(r, t) + 0"(r) (IV.4) 

(/)™ is the gravitational potential field, so that g = —Vcj)"^. We will assume, that (j)'^ does not depend on the time. 

We restrict ourself to systems, where the acceleration of the components relative to each other is small compared to 
the acceleration of the mixture's center of mass. This implies that the kinetic energy of the components relative motion 
is small compared to the kinetic energy of the mixture's center of mass motion. This is true, when the relaxation time 
of the relative motion is very small. For the common mixtures this is the case. Thus, the specific kinetic energy is 

r™(r, i) = lv2(r,t) (IV.5) 



From momentum conservation we obtain: 



dt dxa 

For the internal energy we get 



^ ^ + pv-g (IV.6) 



= -V-Jg - TTa/jV^Q -p V-V + V^ (IV.7) 
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where v^q = dYfj/dxa and where 

Jq = Je-pve™-pv-n-v (IV.8) 

is the total heat flux. 

We write the Gibbs relation Eq. pi.lSp in the form 



i—1 

Using previous equations and performing algebraic transformations we obtain 

p ^ - - V-- (j, - 5] jr) + J.- V- - 5] Ji" • - - TT^^ (IV. 10) 

k=l fc=l 

V. 2D ISOTROPY OF THE SURFACE. 

Even though the fluid does not have any preferred direction microscopically, we cannot say that it has a 3- 
dimensional isotropy everywhere, since there are the mesoscopic directions of the density gradient. The two-phase 
equilibrium state is not 3-dimensionally isotropic. 

A special care should be taken to determine the normal direction to the surface. With the help of the equilibrium 
analysis one can obtain the equilibrium densities distributions in the interfacial region. It is possible therefore to 
determine the equidensity surfaces, i.e. mathematical surfaces, where either density is constant, and which are normal 
to the corresponding density gradient. One may in principal use the gradients of either of the densities to define 
a direction normal to the surface. For the mixture we find it more convenient, however, to define a normal using 
the tension field VQ7Q^(r). We call the surfaces which are everywhere normal to this vector field the equitensional 
surfaces. The thickness of the interfacial region will be assumed to be much smaller then the radii of curvature of 
these equitensional surfaces. Given this assumption the tension vector field in good approximation does not change 
it's direction through the interface. Thus, it is possible to speak about the normal vector n on the surface, which is 
parallel to the tension vector in this region. 

This allows us to speak about the symmetry of the surface. If the surface curvature is the same in both directions, 
parallel to the surface, a small surrounding of the normal through the interfacial region is invariant for any rotations 
around and refiections with respect to this normal. Thus we can say, that such a system has a local 2-dimensional 
isotropy. We shall refer to such a property of the interfacial region as the 2-dimensional isotropy of the surface. If the 
two radii of curvature differ, surface is not 2-dimensionally isotropic any more. For a surface which is thin compared 
to the radii of curvature one can, in a good approximation, consider it to be 2-dimensionally isotropic. We assume 
this to be the case for the systems we will consider. 

If the system has 3-dimensional(3D) isotropy, then coupling occurs only between forces and fluxes of the same 3D 
tensorial character. For an interfacial region, which is 2-dimensionally(2D) isotropic, coupling occurs only between 
forces and fluxes of the same 2D tensorial character. Thus, phenomenological coefficients must remain unchanged 
under rotations and reflections with respect to the direction normal to the surface. Below we show how one can 
extract 2D-isotropic quantities from 3D scalars, vectors and tensors. 

We shall use the special notation for the tensorial quantities of different order and different behavior in this section. 
Any tensorial quantity is denoted as Q^'*''^. Here d indicates the dimensionality of the space, in which the quantity is 
being considered, and can be either 3 or 2 here, r indicates the rank of the tensorial quantity, and can be s for scalar, 
V for vectorial or t for tensorial quantities. We refer to Appendix [B] for the details. 

Gonsider the entropy production, which has a form 

C7s = S'^^'''' R'-^''^ + U^^^^-W^^") + T^^*) : /7(^*) (V.l) 

To be able to use the 2D Curie principle one may proceed along the steps, explained in p^ . To clarify this we 
shall write this expression as a combination of independent 2D scalars, vectors and tensors. The details are given in 
Appendix [B], here we will give the results. 

One can split the vectorial and tensorial quantities into the normal and parallel components with respect to the 
normal vector n on the surface. We use the subscripts ± and || for this quantities. Because of 2D-isotropy of the 
surface, these quantities reveal the scalar, vectorial or the tensorial behavior under rotations around and refiections 
with respect to this normal in a 2D space (we refer to Appendix [B] for the details). This will be indicated by 
superscripts 2 r as explained above. 
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Any 3D scalar is also 2D scalar, since it remains invariant under any kind of coordinate transformations. 

_5,(3s) ^ ^(2s) ^ ^(s) (y 2) 

Any 3D vector V^'^3v) can be written as (cf. the notation with Eq. (jB.ip ): 

y(3v) ^ l^^js)^ ^(2v)^ (y 3^ 

Any 3D tensor T^3t) can be written as (cf. the notation with Eq. ()B.4P ): 
Combining these components we obtain for the entropy production 

(V.5) 

where 

(Js, seal = S^^^ + Vi'^ W['^ + t'i'I n['l + i(TrT|[J*') (Tr 77|fJ*)) 

a,._, = ^^P^).W-|[^^)+2T|^).i7(f^) (V.6) 

_ 7^(2t) . jfi2t) 
^ s, tens II II II 11 

where T^^'' = ^(T||^''' + rj^*-*) and we have used the symmetry of the tensor T^^*). The little circle above a 2x2 

tensor like in T indicates the symmetric traceless part of this tensor. The 2D Curie principle tells us that coupling 
occur only between quantities of the same 2D tensorial order. 

VI. THE PHENOMENOLOGICAL EQUATIONS. 
A. The force-flux relations. 

Comparing Eq. (jIV.lOp with the balance equation for the entropy 

P^ = -V.J, + a, (VI.l) 
at 

we conclude, that the entropy flux and the rate of entropy production are given by 

n— 1 

fc=i 



n-l 



a. =J,-V--^jr-V^--^„^v^„ (VI.3) 

k=l 



According to the second law Cs is non-negative. Comparing Eq. (jVI.SP with Eq. (jV.SP and Eq. (jV.6P we can write the 
entropy production for a 2D-isotropic surface as the sum of 2-dimensional scalar, vectorial and tensorial contributions 



n-l 



CTs, seal = Jg. ± i - E J,"l ^ - (V^V^) ^ - ^ (V|| • V|| ) (Tr TTll || ) 



fc=l 



II T ^ ' 

k=l 

(^s, tens — (V|| V|| 



CTs, veet = J?, H '^H T ^ Yj -^Tw ''^|| ^ 2 V# • y TT^ 
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where = 5(V|[Vj_ + Vj_V||). 

The hnear force-flux equations for the scalar force-flux pairs are 

n-l 
k=l 



n-1 



k=l 



n-1 



^ V_lVj_ — RTTq,±±Jq,± — J2 R-^k,±±Jk,\ ^ RTnr,±±T^±± — ^tttt, ±jl 5 (Tr 7r|[ |[ ) 



k=l 



n-1 

ri^W'^lO = RTrq.\\±Jq,± - J2 RTTk,\\±Jk\ ^ RTTn,\\±T^±± - Rtttt , \\\\ 2 i"^^ ''^ \\\\) 

k=l 



(VI.5) 



For the vectorial force-flux pairs they are 

1 V'" 1 

J9.ll = -^99.1111 ^IIT " E -^9fc, 1111^11+ - ^q^,\\#T^* 

fc=l 

1 0" 1 

k=l 

1 lA" 1 

fe=l 

and for the tensorial force-flux pairs they are 

^,111 = L„i(V||V||) (VI.7) 

All the resistivities R and conductivities L are scalars. One can easily invert the resistivity matrix R and write the 
corresponding relations for the fluxes using the conductivities L, and vice versa. 

For flat surfaces it was found [llj that the resistivities are additive in the normal direction to the surface while 
the conductivities are additive in the parallel direction. We therefore consider it convenient to write the force-flux 
relations in the given form. 



B. The phenomenological coefficients. 

The Onsager relations for the phenomenological coefficients are the following. 



Rqk. ±± 


— Rkq.±± 


Lqk^ nil 


— Lkq^WW 




RqTT,±± 


= RTTq,±± 


Lq^, \\# 


= L^q.#\\ 




Rq-n . _L 


R-Kq. _L 


L^k,\\\\ 


— -^fci, nil 




Rik. ±± 


= Rki. ±± 






(VI 



-^?7r,_L|| -^'Tri, ||_L 

-^7r7T,_L|| -^7r7r,||_L 

As the ordinary Onsager relations, these are the consequence of the microscopic time reversal invariance. 

As usual the values of the phenomenological coefficients locally will depend on the local thermodynamic variables. 
These are the local concentrations, p, and the temperature, T. In the gradient theory the density gradients are 
also considered as local thermodynamic variables. In view of this the phenomenological coefficients may also depend 
on the gradients of the densities. The values of the phenomenological coefficients and their functional dependence 
on the thermodynamic variables are not given by the mesoscopic theory. They should be either calculated from 
statistical mechanical considerations, or from experiments, either real or computer. While they are well investigated 
for homogeneous fluids and fluid mixtures, such data are not available for the surface coefficients of fluid mixtures. 
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We shall use the following expression for each of the resistivity coefficients 

R = R' + {R" - R') + « (R" + fi^) ,J^1' (VI.9) 

where R^ and R^^ arc the resistivities for the coexisting homogeneous phases in the equilibrium state. Here (p is the 
order parameter, typically this is just the density p or the molar concentration c. (p^q is the equilibrium profile and 
jVipeqlmax IS the maximum value of the gradient of this profile. The first two terms are just a smooth transition of 
the resistivity from the value in the one phase to the value in the other phase. This is the first natural assumption 
for the resistivity profile. The origin of the third term comes from the assumption of an excess resistivity in the 
interfacial region. Particularly one can observe this fact in the molecular dynamic simulations [l^ . The exact form of 
this term may be debated. It was chosen to model a rise of the resistivity in the interfacial region. The |Viy9p factor 
makes this term significant only in the interfacial region. It is scaled with |V</3e(j|max order to make this factor 
dimensionless and not far from unity. The {R^^ + R^) factor gives the average value of the resistivity of both phases. 
The dimensionless factor a determines the magnitude of this effect. The homogeneous resistivities R^ and R^^ are 
the known functions of the mass fraction and the temperature along the plane of coexistence. 

For the conductivities used in Eq. (jVI.6p and Eq. (|VI.7p one may use expressions analogous to Eq. (jVI.9[) . The 
conductivities along the surface are expected to be additive Thus it is important to use this equation for the 

conductivities and not for the resistivities along the surface. In this respect it is important to note that a may in 
principle be negative as long as the corresponding R and L remain everywhere positive. For the resistivity this would 
describe an interfacial region with a lower resistivity and for the conductivities it would describe an interfacial region 
with a lower conductivity. Below we will only consider positive a's. 



VII. TYPICAL PROFILES FOR THE BINARY MIXTURE. 



In order to illustrate the results, which one can obtain using the above procedure, we have applied it to a special 
case. This requires a number of approximations, connected with the specific mixture and the geometry. We consider 
a flat liquid-vapor interface of the binary mixture of cyclohexane (f* component) and n-hexane (2"*^ component) in 
non-equilibrium stationary conditions. Wc only consider fluxes and gradients in the direction normal to the surface. 
Furthermore we neglect viscous contributions. The force-flux relations then reduce to 

d 1 p T p J m 

dxT ~ ^qq'^q ^qi^i 

d ip T? J C? Tfn 

_ HiqJq - tCiiJi 

The differential equations for the temperature, density and fraction profiles were solved using a numerical method for 
a two point boundary value problem. This was done using a collocation method implemented in the Matlab function 
bvp4c p^. In this procedure we used the equilibrium profiles, found in Q as an initial state. Further details of the 
solution procedure will be given in the following paper. The numerical values of the homogeneous resistivities were 
taken from In all cases the integrated molar content for both components was kept equal to the equilibrium 

value. 

The first aspect we will try to clarify is the influence of the additional resistivity to transport. We consider in 
particular two cases. In the first only a^g, the a-factor for the heat resistivity coefficient Rqq, is unequal to zero and 
in the second only an, the a-factor for the diffusion resistivity coefficient Ra, is unequal to zero. In Figs. [T][3] we plot 
the total molar concentration, the mole fraction and the temperature for the case that only aqq ^ 0. The system 
is brought out of equilibrium by reducing the pressure on the vapor side to 0.95 Peg, where peq is the equilibrium 
pressure. The temperatures on both ends of the box are kept equal to equilibrium temperature Tg^. In Fig. [1] we 
see, that the total molar concentration does not depend on the value of aqq very much. This is different for the mole 
fraction which increases about 2 % on the vapor side when Uqq increases from to 10. The temperature decreases 
due to the evaporation. In all cases the extrapolated temperature in the liquid is higher then the value in the vapor, 
where we extrapolate to the inflection point of the total molar concentration. For Uqq — 10 the minimum of the 
temperature is below both extrapolations. In that case the temperature "jump" in the the extrapolated profiles has 
increased to about 3 °C. 

In Figs, mo we plot the mole fraction and the temperature for the case that only an ^ 0. Further conditions are 
the same as in the previous example. We did not plot the total molar concentration which is not very different from 
the one given in Fig. [TJ The modification of the mole fraction is now more dramatic. In the vapor it decreases up 
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FIG. 1: Molar concentration profile for different aqq at p = 0.95peg and Tg = Ti = T^q. 

to 27 %. The temperature increases for larger values of an. This is related to a decrease of the evaporation. The 
temperature jump in the extrapolated profiles is in all cases not more than 0.5 °C. 

The third example considers the aqq ^ case when the vapor pressure and temperature are kept equal to the 
equilibrium values and the liquid temperature is 5 % higher then the equilibrium value. The total concentration 
profile is again similar to Fig. [1] and not given. Figs. [6][7] give the mole fraction and temperature profiles. The change 
in the mole fraction went up to 12 %. The temperature jump goes up to about 20 °C. 

As fourth example we consider the an ^ case when the vapor pressure and temperature are kept equal to the 
equilibrium values and the liquid temperature is 5 % higher then the equilibrium value. The total concentration 
profile is still similar to Fig. [T] and not given. Figs. [8][9] give the mole fraction and temperature profiles. One can 
notice again the more dramatic behavior of the mole fraction and temperature profiles for big values of an . 

Next we will discuss the influence of the pressure. We choose both aqq and a\\ equal to 1 for this case. The total 
concentration profiles are again similar to Fig. [T] and not given. In Figs. [TmiTTl we give the mole fraction and the 
temperature profiles. The temperature profile goes down for evaporation and up for condensation as expected. The 
mole fraction in the vapor rises a little bit for condensation and decreases for evaporation. 

In the last case we consider the influence of the temperature on the liquid side. The total concentration profile 
is again similar to Fig. [T] and not given. Lowering (raising) Tg gives evaporation (condensation). This lowers (rises) 
the temperature on the vapor side as expected. The mole fraction decreases (rises) about 5 % for condensation 
(evaporation). This is the opposite of what happens in the previous case. 




VIII. DISCUSSION AND CONCLUSIONS. 



In this paper we have built the framework for the non-equilibrium gradient model for the surface. This require 
the following important steps. With the help of the equilibrium gradient model for the mixtures, established in [J 
we were able to extend the thermodynamic description of the interfacial region to the non-equilibrium case. Explicit 
expressions for the Gibbs relation and the pressure tensor were given. We found, that compared to homogeneous 
systems, the Gibbs relation for the interfacial region contains an additional term, proportional to the divergence of 
the tension tensor. This tensor plays an important role in the surface and integration of it's perpendicular component 
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gives the surface tension. Away from the surface this tensor is equal to zero and we have the famihar Gibbs relation. 
For a single component system the Gibbs relation in the interfacial region reduces to the one given by [1]. 

For the non-equilibrium description we use the standard hydrodynamic equations for the so-called " one-fluid" model 
of the fluid. Together with the Gibbs relation and the balance equation for the entropy density, we were then able 
to obtain explicit expressions for the entropy flux and the entropy production not only in the homogeneous phases 
but also in the interfacial region. This identifies the conjugate thermodynamic forces and fluxes in the interfacial 
region. This made it possible to give the general force-flux relations in this region. The explicit form of these 
equations depends on the symmetry of the system. We discuss why one can consider a fluid-fluid interface to be 
two-dimensional isotropic. Due to the Curie principle coupling occurs only between 2-dimensionally isotropic scalars, 
vectors and tensors of the same rank in such a system. The resulting force-flux relations in the interfacial region are 
accordingly simplified and given. Interesting is that the components of the fiuxes normal to the surface are scalar 
and couple therefore to reactions in the interfacial region. This is related to active transport, a phenomenon of great 
importance. We do not explicitly consider reacting system in this paper, but may do this in a future publication. 
Having all these equations the non-equilibrium description of the surface is complete. 

We applied the description to the special case of stationary heat and mass transport through and into a surface. 
Transport along the surface is in many respects like a 2-dimensional analog of flow in a homogeneous medium. The 
presence of the surface has more influence on transports in the normal direction. We flnd that there is for instance a 
strong effect of the surface on the temperature and the concentration profiles. This effect increases when we increase 
the contribution proportional to the square gradient of the order parameter to the resistivity in the interfacial region. 

We conclude that the description we have given, using an extension of the square gradient model, will be a useful tool 
to study many details of the dynamics of evaporation and condensation in multi-component systems. Non-equilibrium 
molecular dynamic simulations of evaporation and condensation could obtain density, mass fraction and temperature 
profiles. We expect the comparison of these profiles with the present model to be very useful. In particular, this will 
give insight in the size and possible density, mass fraction and temperature dependence of the coefficient a. In this 
manner we expect the model to form a bridge between the microscopic description using non-equilibrium molecular 
dynamic simulations and the discrete macroscopic description using the excess densities introduced by Gibbs flS] . In 
a following paper we intend to investigate whether the discrete description satisfies the local equilibrium assumption 
for an arbitrary choice of the dividing surface, a property which was verified for one-component systems [2]. This 
would be a rather remarkable result, given the fact that the continuous description does not obey this property. For 
a sy stematic development of the non-equilibrium thermodynamic description of surfaces this property is essential 



APPENDIX A: ON THE NON-EQUILIBRIUM GIBBS RELATION. 



For the specific Helmholtz energy the ordinary Gibbs relation in equilibrium was found to be 

n-1 

Sris"^ (r) , v"^ (r) , r (r)) = -s'" (r) + ^ S^T (r) - P (r) Sv"' (r) (A. 1) 

1=1 

Spatial Gibbs relations in equilibrium was found to be 

V/™(r) = y Ver(i-) - P (r) Vw™(r) + w™(r) ^2^£^Il (A.2) 

OXa 

2—1 

One may wonder why we use the Gibbs relations for the specific internal energy, not for the specific Helmholtz 
energy, to extend them to non-equilibrium analysis. Following the same procedure, as in Subsec. |II B| . we can extend 
the Gibbs relations for the specific Helmholtz energy to non-equilibrium in the following way 

i—1 



V/™(r, t) = V ^r(r, t) VCr(r, t)-pir, t) Vt;"(r, t) + v^(v, t) ^22£^ElJl (A.4) 

For the ordinary Gibbs relations, like in homogeneous description, there is no preference in the thermodynamic 
potential. Provided Eq. (jll.4[) . equilibrium ordinary Gibbs relations Eq. pi.7p and Eq. (|A.1|) are equivalent. The 
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non-equilibrium relation between these potentials 



u^{r, t)^r{r, i) + ,s™(r, t)T{r, t) 



(A.5) 



makes non-equilibrium ordinary Gibbs relations Eq. (jll.lip and Eq. (IA.3|) also to be equivalent. 

The situation is different for the spatial Gibbs relations, however. Provided Eq. (|TL4|, equilibrium spatial Gibbs 
relations Eq. (jll.8[) and Eq. (|A.2[) are also equivalent. The non-equilibrium relation Eq. (jA.Sp between these potentials 
makes, however, non-equilibrium spatial Gibbs relations Eq. (jII.12p and Eq. (|A.4|1 to be not equivalent. The reason 
for that is that the equilibrium spatial Gibbs relation for the specific Helmholtz energy does not contain the term, 
proportional to VT, since VT = in equilibrium. In non-equilibrium Vr(r, t) ^ and, as one can see from Eq. (jII.12[) 
and Eq. (jA.sp . V/™(r, t) contains such a term. 

We see, that the Gibbs relations for the specific internal energy describe the system more adequately since they 
do not suffer from the unaccounted effect of possible temperature changes. Because of this reason we should use the 
Gibbs relations for the specific internal energy, not the specific Helmholtz energy, to extend them to non-equilibrium. 

APPENDIX B: 2D ISOTROPIC COMPONENTS IN THE 3D TENSORIAL QUANTITIES 

As in Sec. [V] we shall use the special notation for the tensorial quantities of different order and different behavior 
in this section. Any tensorial quantity is denoted as Q^'''-'. Here d indicates the dimensionality of the space, in which 
the quantity is being considered, and can be either 3 or 2 here, r indicates the rank of the tensorial quantity, and can 
be s for scalar, v for vectorial or t for tensorial quantities. For example, Q^^*) indicates the 2-dimensional tensor, 
i.e. the quantity (^jj 522): where qij are numbers, and Q^^^^ indicates the 3-dimensional vector, i.e. the quantity 
(gi, (72, 93), where qi are numbers. Scalars are the numbers irrespectively of the dimensionality of the space, so they 
will be denoted simply by Q'^^\ 

Some quantities reveal the tensorial behavior of a some rank in d-dimensional space only under some specified 
transformations, while in general they don't. In this section we are interested only in rotations around and reflections 
with respect to some constant vector iV^^v) 3-dimensional space. We will denote quantities which reveal the 
tensorial behavior of rank r under these transformations by Q^'^'^ "\ 

We show how in presence of the constant vector A^'^'^) one can split the tensorial quantity Q^'^'' into a combination 
of the tensorial quantities Q'-^'' Without loss of generality we will assume that N^^ — (1, 0, 0). 

From 3D vector one can construct the following quantities, which are linear in V^(-^^); one scalar quantity 



which is perpendicular to the iV^^v). Denoting ^(2^™) ^ (vj^''\ V^^"'^) we can write that V^^^") = (g, y(2v«))^ 
Thus, 



From 3D tensor T^^'-' one can construct the following quantities, which are linear in T^^^^: 2 scalar quantity. 



and one vectorial quantity 




y(3v) _ yisN) ^(3v) _|_ ^(Svjv) _ fyisN) ^ y(2vN) 



(B.l) 



and 



j.{sn) ^ ^(3v)_y(3t),^(3 



v) _ ^(3t) 
— 11 



two vectorial quantities 



rp{3VN) _ ^(3 V) _rp{3t) _ JiSjy J\^(3v) 



and 
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(which are equal, if T^^*' is symmetric); and tensorial quantity 

/O 

ji(3tN) = 2i(3t) _ J.SJV _/V(3 v) _ j.(3viv) ^(3v) _ jy(3v) 2i(3v) _ |g j^i^t) j^iSt) 

Denoting 



-.(2vn) _ ,'Tn(3t) ^(3t)^ 
; — 1-^12 ' -'is ; 



n(2Vjv) 



(3t) ^(3t)N, 
21 ' -'Sl j 



JT(2tN) 



22 23 
^(3t) ^(3t) 
».-'32 -'-33 J 



we can write that 



j,{3vn) _ ji(2vn)^ 2i(3tN) 



^(3 Viv) 



(0, T.('"")) 








Thus! 

2^(3t) _ j'_(sn) jy(3v) _/y(3v) _|_ j.(3v„) jy(3v) ^ _/y(3v) J^i^^^) _j_ 2^(3tjv) 
j^(2vn) 2,(2tiv) I 

Tensor T^^*™) still contains the scalar part 

j.(sn) ^ rpj. Ji(2tjv) _ rpi^i) I 2^(3*) 

2 32 33 

which obeys the relation 



(B.2) 



(s) _ ^(Siv) _|_ rp(sN) 



(B.3) 



Two of these three scalar quantities are linearly independent and one can use any pair. Since we want to reduce all 
the quantities to the form Q^'^^ we will use T^^^"^ and Tj^"' as independent pair. Introducing the traceless tensor 



j-(2tN) = ji(2tN) _ 1 y(siv) jj(2t) 
2 2 



'^(3t) _ 1 ^(siv) 
-'22 2 -'2 



(3t) 



32 



rp{3t) _ 1 rTn(sjv) 

-'33 2 -'2 y 



and 



jT(3tN) ^ jT(3tN) _ 1 jt(sn) jj(3tjv) 
2 2 








where 



{/(2t) 



1 
1 



[/(3tiv) 





[/(2t) 



we can write a 3D tensor as 



J-CSt) _ y^(siv) ^(3v) jY(3v) _^ rp{3vN) ^(3v) _|_ ^(3v) j^(3vn) _|_ 1 jj(sjv) jj(3tN) _|_ y(3tjv) 



(s,y) 



i(2viv) 



jt(2vn) 1 jT(siv) jj{2t) _|_ 'f{2tN) 



(B.4) 



^ Note, that if the product of two tensorial quantities of rank r > is written without • , it means that this is the product, not the internal 
product. 
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APPENDIX C: HELMHOLTZ ENERGY OF A MIXTURE OF IDEAL GASES. 



According to [l8[ the total Helmholtz energy of a homogeneous mixture of ideal gases is 

F.MT.cA\--RTt..i.{-^^^^^) (C.I) 

k=l '^^ ' 

where Vk is the number of moles and Ck the molar density of component k. Furthermore Afe is the thermal de Broglie 
wavelength and w^ is a characteristic sum over the internal degrees of freedom of component k 

Ak{T) = hNA^2^/MkRT , Wk{T) = J2,eM-4/kBT) (C.2) 

where is the molar mass of component k and are the energy levels of the internal degree of freedom of component 
k. If one describes the mixture using molar specific variables the following equivalent expression is more useful 



(C.3) 



where c is the total molar density of the mixture, 

A(T, = hNA^2Ti/M{0RT (C.4) 

is the mixture's thermal de Broglie wavelength, w(r, f ) a characteristic sum over all the internal degrees of freedom 
of the mixture and 

n n— 1 

fc=l fc=l 

is the molar mass of the mixture. The exact expression for w(T, as well as expression for Wfc(T), is determined by 
model approximation for the mixture. 

The specific Helmholtz energy of a mixture of ideal gases then becomes 



(C.6) 



Due to the spirit of a the onc-fiuid approach we have to equate the second term to 0. Thus, 

/.,.,(r,c,?)^-«Ti„(-£g;|L) (a7) 

where 



w(T.O^exp{|;«.,n(l(-||fw.(T))} 



(C.8) 



can be considered as a mixing rule for the w. We note that Eq. (jC.7[) together with Eq. (|C.8|) does not impose any 
assumptions: it is nothing but Eq. (|C.6p written with the help of one-fluid terms. 



APPENDIX D: SYMBOLS LIST 





- contraction sign 




- double contraction sign 


1, n 


- enumeration of all integers from 1 to n 


Tr 


- trace 


f 


- traceless part of a tensor 




- cartesian indices 


± 


- perpendicular direction to the surface 


II 


- parallel direction to the surface 


d 


- partial differential 


d,5 


- differential 


Sap 


- Kroneker symbol 
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V 


[1/m] 


- 


nabla operator 


7a/3 


[ Pa ] 


- 


tension tensor 




[ J/kg ] 


- 


square gradient coutribution 




[ J mVkg^ ] 


- 


square gradient coefficients for molar units 


Ki 


[ J m/kg ] 


- 


>> >> 


kV^ 


[ J/m^ ] 


- 


>> >> 




[ J mVkg^ ] 


- 


square gradient coefficients for volume units 


A 


[ ™ ] 


- 


thermal de Broglie wavelength 




[ J/kg ] 


- 


mass chemical potential of component n 


n 


[ Pa ] 


- 


viscous pressure tensor 




[ Pa ] 


- 


>> >> 


P 


[ kg/m^ ] 


- 


mass density 


(h. 


[ kg/m^ ] 


- 


mass density of component i 


CTs 


[ J/(K s) ] 


— 


entropy production 


0'al3 


[ Pa ] 


- 


total pressure tensor 




[ J/kg ] 


— 


kinetic energy density per unit of mass 




[ J/kg ] 


- 


potential energy density per unit of mass 




[-] 


- 


molar fraction of component i 


Si 


[-] 


- 


mass fraction of component i 




[ J/kg ] 


- 


reduced mass chemical potential of component i 


A 


[ kg / (mol s)^ ] 
[ kg mV (mol sf ] 


- 


van der Waals equation of state coefficient for a mixture 


tti 


- 


van der Waals equation of state coefficient for a pure component i 


B 


[ m'^/mol ] 


- 


van der Waals equation of state coefficient for a mixture 


bi 


[ m' /mol ] 


- 


van der Waals equation of state coefficient for a pure component i 


c 


[ mol/m'^ ] 


- 


molar concentration 


Cm 


[ J/kg ] 


- 


total energy density per unit of mass 


r 


[ J/kg ] 


- 


Helmholtz energy density per unit of mass 


Jo 


[ J/kg ] 
[ m/s^ ] 


- 


homogeneous Helmholtz energy density per unit of mass 


g 


- 


gravitational acceleration 


i,j,k 


[-1 


- 


component number 


Je 


[ J/(m^s) ] 


- 


total energy flux 


■f m 


[ kg/(m2s) ] 
[ J/(m2s) ] 


- 


total mass flux of component k 


J. 


- 


total heat flux 




- 


phenomenological conductivities 


Na 


[ 1/mol ] 


- 


Avogadro's number 


n 




- 


number of components 


P 

Q(dr) 


[ Pa ] 


- 


pressure 


[ ~ ] 


- 


tensorial quantity of rank r in d-dimensional space 


R 


[ J/(K mol) ] 


— 


universal gas constant 


Rah,xy 


[~1 


- 


phenomenological resistivities 


r 


[m] 


- 


position 




[ J/(K kg) ] 


- 


entropy density per unit of mass 


T 


[K] 


- 


temperature 


t 


[s] 


- 


time 


vT 


[ J/kg ] 


- 


internal energy density per unit of mass 


V 


[ m^/mol ] 




molar volume 


v"" 


[ m'Vkg ] 




volume per unit of mass 


V 


[ m/s ] 




barycentric velocity 


Vi 


[ m/s ] 




velocity of component i 


a; 


[m] 




position 
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FIG. 2: Molar fraction profile for different aqq at p = 0.95peg and Tg =T( = T^q. 
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FIG. 3: Temperature profile for different aqq at p = O.QSpeq and Tg = Tt = Teq. 
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FIG. 4: Molar fraction profile for different an at p = 0.95pe 
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FIG. 5: Temperature profile for different an at p = 0.95peg 
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FIG. 6: Molar fraction profile for different aqq at p = peq, Tg = Teq and Te = l.OSTe, 
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FIG. 7: Temperature profile for different Oqq at p = peq, Tg = Teq and Tt = 1.05Teq 
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FIG. 8: Molar fraction profile for different an at p = peq, Tg = Teq and T( = l.OSTeq 
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FIG. 9: Temperature profile for different an at p = peq, Tg = Teq and Tt = 1.05Teq 
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FIG. 10: Molar fraction profile for different pa.tTg=Te = Te, 
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FIG. 11: Temperature profile for different p at Tg = Te = Teq 
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FIG. 12: Molar fraction profile for different Tg, at Tg = Teq and p = Peq 
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FIG. 13: Temperature profile for different Ti at Tg = Teq and p = peq 



